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ABSTRACT 

Techniques for extractions boundary layer parameters from measurements 
of a short-pulse (w 1 fxs) CO 2 Doppler Lidar (A = 10.6 fim) are described. 

The Lidar is operated by the National Oceanic and Atmospheric Administration 
(NOAA) Wave Propagation Laboratory (WPL). The measurements are those col- 
lected during the First International Satellites Land surface Climatology Project 
(ISLSCP) Field Experiment (FIFE). The radial velocity measurements have a 
range resolution of 150 m. With a pulse repetition rate of 20 Hz, it is possible to 
perform scannings in two perpendicular vertical planes (x-z and y-z) in « 72 s. 

By continuously operating the Lidar for about an hour, one can extract stable 
statistics of the radial velocities. Assuming that the turbulence is horizontally ho- 
mogeneous, we have estimated the mean wind, its standard deviations, and the 
momentum fluxes. From the vertically pointing beam we have estimated the first, 
second, and, third moments of the vertical velocity. Spectral analysis of the ra- 
dial velocities is also performed from which, by examining the amplitude of the 
power spectrum at the inertial range, we have deduced the kinetic energy dissipa- 
tion. Finally, using the statistical form of the Navier-Stokes equations, the surface 
heat flux is derived as the residual balance between the vertical gradient of the 
third moment of the vertical velocity and the kinetic energy dissipation. With the 
exception of the vertically pointing beam, an individual radial velocity estimate 
is accurate only to ±0.7 ms'* 1 . Combining many measurements would normally 
reduce the error provided that, it is un-biased and un-correlated. The nature of 
some of the algorithms however, is such that, biased and correlated errors may be 
generated even though the “raw” measurements are not. We have developed data 
processing procedures that eliminate bias and minimize error correlation. Once 
bias and error correlations are accounted for, the large sample size is shown to 
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reduce the errors substantially. We show for instance that, a single momentum 
flux estimate has an accuracy of ±2 m 2 s -2 but, when combined with other mea- 
surements the error can be reduced to ±0.10 m 2 s -2 . The principal features of the 
derived turbulence statistics for two case studies (11 July 1987 16:11-17:10 UTC 
and 17:29-18:10 UTC) are as follows: 

(a) The mean surface wind is from the south and has a speed of ~ 15 ms -1 . 
There is a southerly jet at 1 km above the surface. The wind in the direction nor- 
mal to the surface wind is becoming more westerly with height. 

(b) The derived momentum fluxes in the mixed layer agree with the unfiltered 
airplane measurements. They are of the order of « 0.5 m 2 s -2 near the surface 
and are retarding the southerly wind. At the stable layer, the momentum fluxes 
are counter-gradient; they remain negative even though the southerlies are di- 
minishing with height; and furthermore they are concentrated in thin layers. By 
spectrally analysing the vertical beam, we have identified the Brunt Vaisala fre- 
quency as the dominant frequency in the stable layer. The available evidence sug- 
gests that the large counter-gradient fluxes are due to a critical layer singularity. 
Under these conditions gravity waves flux divergence may be large. Current Gen- 
eral Circulations Models (GCM’s) generally do not parameterize counter-gradient 
momentum transport. If the phenomenon is prevalent, our study suggests that in- 
clusion of such an effect may be important. 

(c) The spectrum of the wind near the surface has a —5/3 power law at the iner- 
tial range. The heat flux estimate is 100 ± 20 Wm' 2 in apparent agreement with 
the surface stations estimate but, systematically larger than the filtered airplane 
data. 
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1. Introduction and Goals 

The purpose of this paper is two fold. First we wish to describe in some 
detail an analysis technique of a single Doppler Lidar (for Light Detection and 
Range) data which allows the estimation of the following turbulence parameters: 
a) mean wind (u, v, w) where ( ) is a symbol for some kind of average (to be de- 
fined in section 4) and (u, v, w ) are the Cartesian velocities in the (x, y, z) direc- 
tions respectively; b) the variance (u' 2 ,u' 2 ,u/ 2 ), where ( )* is a symbol for the 
deviations from the mean; c) the part of the covariance associated with the ver- 
tical fluxes of horizontal momentum ( u'w ', v'w ') and d) the third moment of the 
vertical velocity, w' 3 , the kinetic energy dissipation e and the surface heat flux. 

A second objective of our research is to study the behaviour of the convective 
Planetary Boundary Layer (PBL) and it’s overlying stable region under the con- 
dition of strong surface wind, moderate surface heat flux and a low level jet. We 
shall sometimes refer to the PBL as ABL (Atmospheric Boundary Layer). 

The scale of motions investigated in our study are in the range of 
150 m - 12 km, too small to be resolved by current and near future General Cir- 
culation Models (GCM’s); yet generally recognised to have an important influ- 
ence on the larger scales. The proper “parameterisations” of these processes re- 
mains an important unsolved issue despite several decades of research (Wyngaard, 
1983). One reason for the apparent slow progress is the requirements for acquir- 
ing large samples in a relatively short time so that, on the one hand the statistics 
are stable but, on the other hand the sampling time must be sufficiently short, 
so that evolution can be ignored (Wyngaard, 1983). As a result “too many mod- 
els are chasing too little data...”. Due to an ability to measure simultaneously 
a ray of points (actually resolution volumes) remote sensing by Doppler Lidars 
and Radars could contribute significantly to alleviate the sampling problems . In 
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our case for instance, a ray (containing fifty to a hundred data points) could be 
sensed in 1/20 s. Such sampling rate is an order of magnitude faster than in situ 
observations (e.g by aircraft or meteorological towers). 

Another unresolved issue in atmospheric turbulence is the similarities and 
differences between stratified and convective turbulence (Lilly, 1982) as well as 
the interaction between the convection in the mixed layer and the waves in the 
stable region above (Kuetner et al., 1987). The Lidar data that we have collected 
include observations in both the mixed layer and parts of the overlying stable re- 
gion. Thus, it is possible to compare theories with observations. 

The data base for this study are measurements taken during the First Inter- 
national Satellites Land Surface Climatology Project (ISLSCP) Field Experiment 
(FIFE) in Manhattan, KS. Since, the experimental technique, the data reduction 
and the hardware are still under development (Eberhard et al., 1989) we confine 
ourselves here to two case studies (11 July 1987, 16:11-17:10 UTC and 17:29-18:20 
UTC). The Central Standard Time (CST) is lagging five hours behind the Coor- 
dinated Universal Time (UTC). 

2. System Characteristics 

For the purpose of this presentation we assume that the Doppler Lidar is 
simply a “black box” which measures the frequency shift between a transmitted 
(coherent) pulse and the backscattered power. The frequency difference is a mea- 
sure of the velocity component toward or away from the Lidar (i.e., the radial ve- 
locity). The frequency u ; of the transmitted pulse is given by u> = c2n / \ where 
A is the wave length of the emitted radiation (10.6 fim) and c the speed of light. 
The pulse duration AT is 1/^s which means that at any given time the recorded 
backscattered power is actually a superposition of returns from several ranges re- 
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fleeting different parts of the pulse. The range resolution A R (i.e., the length of 
the line interval which corresponds to backscattered power recorded at the same 
time) is then given by 2A R = cAT. This means in our case that A R = 150 m. 
The factor 2 takes into account that the pulse goes out and is reflected back. 

The minimum range below which velocity estimates are not reliable is 600 m. 
The maximum range above which the signal/noise ratio is too low to permit ra- 
dial wind estimates depends mostly on aerosols content and absolute humidity 
(Eberhard et al., 1989). For these particular case studies it appears to be 6 km in 
the horizontal and 2 km in the vertical. Since the Lidar scanner can move ver- 
tically and horizontally, one can combine “forward” and “backward” shots to 
obtain at the surface velocity measurements along a 12 km line with a “gap” of 
seven range gates in the centre. As one moves upward the lines becomes progres- 
sively shorter since no useful returns are available above 2 km. Finally it should 
be mentioned that the pulse repetition rate is 20 Hz (i.e., 20 pulses are being 
fired in one second). With the exception of the vertical direction, three pulses 
are combined to get a more stable (less noisy) velocity estimate (±0.70 ms -1 ). 

In the vertical direction we use twelve shots to get a more accurate estimate 
(±0.35 ms -1 ) of the vertical velocity . 

Combining the above information it is concluded that, a radial can be sensed 
in 0.15 s. Since a single horizontal beam contains 37 useful velocity measure- 
ments the over all sampling rate is several orders of magnitude faster than typical 
sampling rate of an airplane. Further details of the Lidar characteristics may be 
found in Eberhard et al. (1989). 

3. Experimental Design 

The experimental design is a compromise between two conflicting require- 
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ments. First one a wishes to cover every resolution volume in the boundary layer 
(i.e., to observe a contiguous sequence of volumes, each approximately 150 m 3 in 
size. Second, it is deemed desirable to scan the entire boundary layer fast enough 
so that, the assumption of a statistical steady state would be valid. Unfortu- 
nately, taking into account the system characteristics it is not possible to satisfy 
the above requirements simultaneously. Instead of covering the entire volume we 
have decided to scan in two mutually orthogonal planes. For horizontally homo- 
geneous turbulence the choice of the two orthogonal planes could be arbitrary; 
knowledge of the statistics in one coordinate system implies knowledge in any 
other coordinate system (Batchelor, 1956). For convenience however, one plane 
say (x-z) is formed by the vertical and the line in the direction of the mean sur- 
face wind. The other plane (y-z) is the vertical and the direction orthogonal to 
the mean surface wind. This type of scanning precluded an estimation of u f v * , 

To accomplish the scanning described above an a’priori knowledge of the mean 
surface wind need to be assumed and this can be obtained from a single Velocity 
Azimuth Display (VAD) scanning (Kropfli, 1986). 

The actual implementation of the scan (called also RHI for Range Height In- 
dicator) consists of first positioning the Lidar (at zero elevation) in the direction 
of the mean surface wind. Keeping the horizontal direction fixed, increments in 
elevation 0 are made until an elevation 0i is achieved. Both 0i and A0 (the in- 
crement) are controllable. In our case we have taken A0 ^ 1° and 0i =42°. 
From the position 0 j one can move the Lidar as quickly as possible to a vertical 
scanning, which is ±2° away from the zenith .This is followed by a quick descent 
to 0! . From 0! the Lidar moves down to zero elevation by decrements of A0. 
This completes one RHI scan. Next, one points the Lidar in a direction orthogo- 
nal to the first RHI plane and performs another RHI scan. At the rate of 20 Hz it 


8 


is possible to complete the above described two RHI’s in ~ 72 s. Whenever pos- 
sible, such a scanning mode is performed uninterruptedly, for about an hour so 
that, stable statistics may be obtained. 

Figs. 1 and 2 are coloured graphical displays of the radial velocities during 
one of our RHI scannings. Red colours are Doppler velocities pointing away from 
the Lidar (red shift) and blue colours are when the targets are moving toward 
the Lidar. The full colour scale is indicated in the figures. Fig. 1 is for an RHI 
plane with the x-axis pointing in the direction of the surface wind and Fig 2 is in 
a plane normal to the first plane. In Fig. 1 the velocities are positive (red shift) 
when the Lidar is 5° away from the north and are negative (blue shift) when the 
shots are in the direction of 185°. This means that mean surface wind is essen- 
tially from the south. A low level jet at an height approximately one kilometre 
above the ground is also prominent. In Fig. 2 the eddy structure is apparent as 
well as a tendency of the wind to become more westerly with height. This latter 
property is revealed most prominently at the jet level where where red shifts are 
apparent in the positive y axis (azimuth 95°); blue shifts are recognized in the 
opposite direction. 

4. Analysis Technique 

a. Mean wind 

Let a denote a radial velocity (positive if it moves away from the Lidar ) and 
let u, v, w denote Cartesian velocities. Let x, y, z be Cartesian right handed coor- 
dinates with the positive x being in the direction of the mean surface wind. Let 0 
be an elevation angle. From geometry we find, 
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a = u cos 0 ± w sin 0 ; LIDAR pointing in the direction of the mean wind . (1) 
a = — u cos 0 + w sin 0 ; LIDAR pointing opposite to the mean wind . (2) 

a = v cos 0 + w sin 0 ; LIDAR pointing 90° normal to the mean wind . (3) 

a = —v cos 0 ± w sin 0 ; LIDAR pointing 270° normal to the mean wind . (4) 

Assuming momentarily that u, v and w are constants at fixed height and el- 
evation angle, it is apparent that, u and w can be found in the x-z plane while, 
v and w may be estimated in the y-z plane. In practice of course, u, v, and w 
are not constants. In that case let t be an index for observations which are con- 
strained to be sufficiently close (say ±50 m) to a certain height zo (e.g., z 0 = 

500 m or zq = 600 m etc. ). Let u,v,w be the ensemble average of the veloci- 
ties in the above described “bin” . The average is both temporal and spatial. We 
may then rewrite (1), (2), (3), (4) as: 


a, = ±u cos 0; + w sin 0, + fluctuations ; i = 1, 2, ...., N\ ; along the wind . (5) 


a{ = ±u cos 0^ + re sin 0, + fluctuations ; i = 1, 2, ...., A ^2 ; across the wind. (6) 
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Here N j and iV 2 axe the number of observations in the along and across the wind 
“bins” respectively, corresponding to a certain height z 0 . The (+) comes when 
the Lidar points in the positive x or positive y directions respectively. The nega- 
tive sign is for negative directions. 

Relations (5) define Ni equations with two unknown (u,w) which can be 
solved by conventional “least square” techniques. Because only the l.h.s of (5) 
is known the fluctuations are omitted. For horizontally homogeneous turbulence 
and non-biased observations this approach may be formally justified i.e., 


N\ —* oo =*> ( u , w ) — > (ensemble average of u and w) . 

Likewise, under similar assumptions (u, w ) can be found in the y-z plane. 

In practice the above least square method is slightly modified. The u profile 
is estimated from data where most of the beams have a low elevation angle. De- 
pending on the application two sets of w are generated, one solely from the ver- 
tically pointing beams (±2° away from the zenith) and, the other from inclined 
beams with elevation angles in the interval 30° < © < 42°. 

It should also be mentioned that, at low heights most of the observations cor- 
respond to low elevation angles. In that case the estimate of w is not reliable. 
With the closest range being 600 m we have found that w cannot be computed 
reliably below 400 m . For such low heights a different method for estimating 
the ensemble average vertical velocity has been devised. For low elevations (say 
0 < 10°) and assuming as a rule of thumb that, u u)andu w , we write for a 
range R and an elevation 0. 


a, = ±ucos 0,+ 


du 

dx 


R cos 2 0 , + fluctuations 


1,2, ....,iVi ; along wind. (7 — a) 
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a, 


dt) 

±v cos 0, -f — R cos 2 0,+ fluctuations ; i = 1, 2, 

oy 


Ni ; across wind . (7 — 6) 


Relation (7) and (8) may be viewed as modifications of (5) and (6), where 
the mean gradients are being represented explicitly rather than absorbed in the 
fluctuations. The mean vertical velocity for low elevations is then obtained from 
vertically integrating the Boussinesq’s form of the mass continuity equation (e.g., 
Ogura and Phillips, 1962) 


w — w 0 — f (du/dx + dv/dy) . (8) 

Jo 

Here wq is the mean vertical velocity at the surface. In this study we assume 
that, wq = 0 . 

b. Second moments 

To estimate the turbulence statistics we decompose the Cartesian velocities 
in the usual way as, 


u = u + u' ; v = v -f v' ; w = w + w' ; a = a + a' . 

We substitute the above in (1), (2), (3), (4) and, define the mean radial velocity 
as, 


a = ±ucos0, + wsinQi along wind ; a = ±ucos0,- + wsinQi ; across wind . (9) 
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Using straight-forward algebra and the assumption of horizontal homogeneity we 
a bin corresponding to an height z 0 that, 

Along the wind 

(o' ), = ±u'w'sin2Qi + w' 2 sin 2 ®i + u' 3 cos 2 Q{ ;i = l,2...,iVi , 

+ fluctuations. 


Along the wind 

(o' )i = ±v'w'sin2®i + w' 2 sin 2 ®i + v' 2 cos 2 ®i ; i = 1,2 , 

+ fluctuations. 

From the above it can be deduced that, the second moments u' 2 , v' 2 , w' 2 , u'w 1 , 
v'w' can be obtained by techniques similar to the one deployed for the calcula- 
tions of the mean wind i.e., a least square solution of N\ ( or iV 2 ) equations. 

Similar to the mean wind estimations (considered in the previous section) 
the reliability of the results depend on favourable geometries. Consequently, for 
each moment calculations different data stratification are used. The variance of 
the vertical velocity has 1 /sin 2 Q dependency. Therefore, at low heights (where 
most of the data comes from too low elevation angles) the w' 2 estimate is sus- 
pect. Likewise, the horizontal momentum fluxes u'w' and v'w' have a l/sm20 
dependency. Thus, the bulk of the observations used for estimating the momen- 
tum fluxes must come from intermediate elevations say, 30° < 0 < 60° . The 
horizontal velocities variances varies ~ 1 /cos 2 ® and consequently most of the 
data must be at low elevations. Taking into account the above constraints as well 


( 11 ) 


( 10 ) 
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as the range limitations in the horizontal and vertical directions, we find that the 
momentum fluxes cannot be estimated reliably at heights lower than 400 m. The 
vertical velocity variance can be reliably estimated only from the beams which are 
less than 2° away from the zenith (i.e., the analysis is limited to heights greater 
than 600 m). The variance of the horizontal velocities on the other hand is esti- 
mated from the ground up to an height close to 2 km. 

To get the vertical velocity variance at levels below 600 m one needs to inte- 
grate the continuity equation from the bottom up, setting zero vertical velocities 
as the lower boundary conditions. It is important to remember that only the ra- 
dial winds are given. Thus, u and v are not measured at the same points. We 
must therefore, consider the statistical form of the continuity equation. The de- 
terministic equation for the fluctuations is given viz., 

w'(z ) = — J ( du'/dx + dv'/dy^jdz'. 

Therefore, the vertical velocity variance is, 




// 

0 0 


du^z') du f (z ") dv^z 1 ) dv*{z n ) dv f (z ff ) 


dx 


dx 


dy dy 


dx 


dy 


dz’dz" . 


Assuming the turbulence to be horizontally isotropic we obtain, 


2 // 


du , (z t ) du\z n ) 
dx dx 


o o L 


dz'dz n ; along the mean wind. (12) 


z z 

= 2 // 


o o 


dv'(z’) dv\z^) 
dy dy 


dz'dz" ; across the mean wind. (13) 
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In practice, instead of calculating the double integrals in (12) and (13), a 
computationally simpler approximate procedure is adopted. For low elevation an- 
gles, 


1 d du 

RBR RV ^d^'’ :miVe ^ W - 


(14) 


Here v e is the meridional velocity (i.e., the tangential velocity in the cylindrical 
coordinates R, 0) and v R is the radial velocity. Define now w' and v' Q as the so- 
lution (with mean quantities subtracted) of the 2-D continuity equation with zero 
vertical velocity at the bottom namely, 


1 d n , Id, 

R d. R RVr + R 30 ”® - 


(15) 


Let (as before) i be an index for a “bin” i.e., all the observations with height 
z, t above the ground and satisfying the condition, 


z, k — A z < Zi < Zi k + A z . 
From (12)-(16) we conclude that, 


(16) 


< = 2 /n ( 17 > 

u = l 

At present we are not yet able to obtain definitive vertical velocity variance 
from the above vertical integration method. Therefore, they will not be presented 


here. 
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c. The vertical velocity balance and the surface heat flux. 

The balance of vertical velocity fluctuations for horizontally homogeneous 
turbulence may be shown to be (Wyngaard and Cote, 1971) as, 


E<5 5 ’’ i*i=* 


(18) 


Here p f is the pressure fluctuations, p Q is the air density at the surface, 
e is the kinetic energy dissipation, g is the gravity acceleration, d 0 is the 
potential temperature at the surface and $* is the potential temperature 
perturbations. Other symbols have been defined elsewhere. Since w *0* is 
the potential temperature — vertical velocity covariance it follows that, p o C p w i 0 f 
is the heat flux with C p the heat capacity of air at constant pressure. 

Wyngaard and Cote (ibid.) find that under unstable conditions (i.e., positive 
heat flux) the pressure term at the surface layer (say the lowest 50 m of the PBL) 
is small compared to the other three terms in (18). Therefore, at the surface layer 
an estimation of the kinetic energy dissipation and the vertical gradient of the 
third moment of the vertical velocity would enable the calculations of the surface 
heat flux as a residual. It should be emphasized that the vertical velocity pressure 
correlation is (according to Wyngaard and Cote) an important term in the kinetic 
energy balance. On the the other hand at the surface layer, the correlation be- 
tween the vertical velocity and the vertical pressure gradient is small compared to 
the other terms in (18). 

Now, a traditional way to estimate the kinetic energy dissipation is by exam- 
ining the line spectra of the longitudinal velocity correlation. In the inertial range 
the expected relation is (Batchelor, 1956), 
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/(«) = ae 2/3 « _5/3 . (19) 

Where k is the wave number, a a universal constant (= 0.5) and /(k) the 
Fourier transform of the longitudinal velocity correlation (i.e., the transform of 
u'(x)u'(x + r) or v'(y)v'(y + r) . 

As has already been mentioned in the experimental design section, direct 
measurements of vertical velocities (from a vertically pointing beam) are available 
for heights greater than 600 m. This data is used to calculate the third moment 
of the vertical velocity w' 3 . For heights lower than 600 m one needs to integrate 
the mass continuity equation in ways similar to what have been used in the esti- 
mations of the second moments. The analogues of (12) and (13) turn out to be, 


“’" w = 2 /// 


du'(z') du'(z ") du'(z'") 


dx 


dx 


dx 


ooo 


dz'dz"dz m ; along wind . (20) 


= 2 /// 

0 0 0 

In practice, one does not need to evaluate explicitly the triple integrals in 
(20) and (21). Instead, the two dimensional continuity equation (14) with the 
provisions stipulated by (15) and (16) is vertically integrated. If the turbulence 
is horizontally isotropic we obtain, 


Ov'(z') dv'{z") dv'(z"') 


dy dy 


dy 


dz'dz"dz"' ; across wind . (21) 


iv 


<V; 

r , = 2 / lV i E 

ifc = l 




( 22 ) 
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Definitive results from the above integration method are not yet available. 
Consequently, in this paper we limit our discussions to third moments data at 
heights greater than 600 m. 

5. Error Analysis 

a. General considerations 

In the discussions below we confine ourselves to the effects observational er- 
rors have on the estimates of the mean wind and its various moments. Other er- 
ror sources include inter alia, insufficiently large number of independent obser- 
vations (Wyngaard, 1983) and inappropriate modelling, e.g., the turbulence may 
not be horizontally homogeneous. For an analysis of the vertically pointing beam 
data there is an additional source of error namely, a Taylor-like hypothesis which 
presumes that time-averaged statistics are the same as line-averaged statistics. 
This latter assumption is clearly dependent on the scales of motion and the aver- 
aging time; “large” eddies may move sufficiently slow and thus may not be ade- 
quately sampled. Our sample size appears to be adequate. We find for instance 
that, experimentations with subsets of the sample (e.g. omitting every other vol- 
ume scan or every other range) have resulted in statistics which are quite close to 
those obtained from the complete sample (i0.1m 2 s -2 for the second moments). 

In discussing the sensitivity of the results to observational errors a distinc- 
tion must be made between biased and unbiased errors as well as correlated vs. 
uncorrelated errors. Let /?* stands for some observation [e.g., square of the radial 
velocity fluctuations in (10) and (11) or the radial velocities in (5) and (6)]. Let 
Sdi be the observational error. Let { } be a symbol for the expected value. The 
errors are considered unbiased if, 
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(W = 0 . (23) 

The errors are uncorrelated (or white) if all the off-diagonal elements of the error 
covariance matrix are zero i.e, 


mtfij) = CiA i • (24) 

Where <5,j is the Kronecker delta and c lt are constants. 

We shall also presume lack of correlation between the exact value /?, and the error 

t>Pi i-e., 


OWi) = 0 (25) 

Now, for unbiased white noise it is known that, as long as the overdeter- 
mined systems of the type (5), (6), (10) or (11) are well conditioned ( i.e., small 
changes in the r.h.s lead to small changes in the solution) the estimates become 
more accurate as the number of independent observations increase. Furthermore, 
the general il rule of thumb” is that, if a(Pi) is the standard deviation of the er- 
ror of an individual observation, the error estimate of the solution is proportional 
to a{Pi)/\/N where N is the number of independent observations. By contrast, 
when bias is present increasing the number of observations does not in general, 
reduce the errors. The above remarks are intuitively obvious when the estimates 
are of the form aj3 t with £)a t - = 1- However, their range of applicability is 
much broader (Jazwinsky, 1970). 

The analysis becomes more involved when the errors are correlated. In fact, 
the simple least-square approach considered in this work is no longer optimal 
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(Jazwinsky, 1970). In the simple approach one tries to minimize a functional J 
defined as, 



( 26 ) 


With a t] some matrix and x } the vector of the unknowns. In the more elaborate 
approach one minimizes the functional J defined as 


J = ^2 ^2 ^2 XI Kim(ais x s ~ 0i){a mt x t ~ 0m) ■ (27) 

mist 

Here Ki m is the covariance matrix of the observational errors i.e., 


K im = (SfcSfr) . (28) 

Of course, in the special case where K{ m = c8 im (i.e., non correlated uniform 
observational errors) (27) is reduced to (26). 

In what follows, we limit our analysis to non-correlated errors. As is shown 
further below the computational procedures that we have developed reduce con- 
siderably the error correlations. 

While not entirely justified, we assume that, the “raw” radial velocity ob- 
servations are not biased. An important source of bias for CO 2 Doppler lidars 
is a systematic deviations from the reference frequency. To partially overcome 
this problem we have eliminated all observations where the pulse monitor signal is 
weak or displaying a tendency toward positive or negative bias. We have also in- 
spected the data manually. After some experimentations we have decided (subjec- 
tively) to, remove all data points where the absolute value of the radial velocity is 
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larger than 25 ms -1 . For each beam we have computed an average and then, re- 
moved all data points whose deviations from the average exceed in absolute value 
7 ms -1 . Virtually identical statistics are obtained if the cutoff value is 5 ms -1 . 
Another source of bias is pointing accuracy. Eberhard et al. (1989) discuss the 
(subjective) correction procedure for this type of errors. 

b. Mean wind 

According to Eberhard et al. (ibid.) the random error associated with an 
individual radial velocity estimate (using 3 shots) is ±0.7 ms -1 . Typically more 
than 2000 observations are used to estimate u and v. Therefore, the uncertainty 
of the estimate of the mean horizontal wind due to uncorrelated non-biased obser- 
vational error is quite small (±1.5 cm s -1 ). 

Recall that in one method of computations the mean vertical velocity above 
600 m is estimated from the vertically pointing beams. Also note (section 2) that 
12 shots are fired in the vertical direction, thus the error of individual radial ve- 
locity estimate is ±0.35 ms -1 . The sample size in the vertical is is also quite 
large. If, all beams which are less than 2° away from the zenith are included we 
have N ~ 1800. Therefore, the uncertainty of w due to non-correlated errors 
is ±1 cms -1 . In another method of computations the w is estimated from data 
where the elevation angle is in the interval 30 < 0 < 42. With a smaller sample 
size (N=600) less favourable geometry (0 = 36) and an uncertainty of individual 
measurements ±0.7 ms -1 , this latter method has ±5 cm/s. On the other hand, 
this latter estimate may be more representative than the former. 

c. Momentum flux 


Let’s 0*. be a single elevation which belongs to a “bin” (i.e., a vertical layer) 
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k. Let m be an index for a volume scan and let v' r be the radial velocity devia- 
tions from the mean. We can rewrite (10) or (11) as, 


( v' r ) fc m .j = u' 2 cos 2 0*, + w' 2 sin 2 0*, + u'w' sin 20*, . ;for positive x . (29) 


And, 


( v 'r ) k m-2 = u ' 7 cos2 d* w>i S ^ n2 — u'®' sin 20*. . ; for negative x . (30) 

Here (; 1) is an index for the forward beam and (;2) an index for the backward 
beam. 

Subtraction of (30) from (29) yields, 


u'w' k . m (individual measurement) — 



«\ m ;2]/( sin20 ‘.)- (3D 


Summing over all volume scans and all elevations within a single “bin” k one ob- 
tains, 


u'u.'Jall measurements) = - (■v ) i , m;2 ]/ ( sm20 t.j . 


(32) 


An important point to emphasis is that, an individual v' r (or v' 2 ) may be 
biased. However, due to subtraction and the homogeneity of the errors (there is 
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only one instrument) the flux error Su'w 1 is not. Furthermore, there is no rea- 
son to believe that, the off-diagonal terms of the error covariance matrix are not 
zero. In addition according to (25) little or no correlation is supposed to exist be- 
tween a radial velocity error 8v r and the “true” radial velocity v r . Under these 
conditions we can estimate the mean and standard deviation of 8u'w' using con- 
ventional error analysis. 

We start by noting that, 


8 (v?) = 2v' r 8v' r + (8v' r ) 2 . (33) 

From the error homogeneity for two different observations with indices £,m; 1 and 
k,nr, 2 [see (32)] we can deduce (formally) that, an individual flux estimate is un- 
biased because, 

($( v r — £( v r )fc,m; 2^ = 0. (34) 

Utilizing (33), (34) and the various non-correlation properties discussed above we 
obtain, 

([^(^r )fc,rn;l - <5(^r )t- i m; 2 ] 2 ) 7 = 2 3/ 2 a(v' r )a(8v' r ) . (35) 

Where the reader is reminded that a is the standard deviation. 

From (32)-(35) and the various error properties it is not difficult to deduce 

that, 


(Su'w 1 ) = 0 . 


(36) 



cr(Su'w') = 2 1 / 2 N l t 2 1 sin(20) cr(v' r ) cr(6v' r ) . 


(37) 


Here 0 is some mean elevation angle for the “bin” k. Since the flux retrieval is 
limited to elevation angles between 30° to 42° we take 0 = 36°. We assume 
(based on the horizontal homogeneity hypothesis) that, 


cr 




(38) 


We find cr(uj.) « ±2 ms 1 and from Eberhard et al. (ibid.) a(Sv ' r ) = ±0.7 ms *. 
Typical values for N are ~ 600. Over all then, a(Su'w l ) = ±0.1 m 2 s -2 

d. Variance and third moments 

As has already been indicated straight-forward squaring of the radial veloc- 
ity induces bias and error correlation. Inspection of (10) and (11) reveals that 
the variance calculations, u' 2 , v' 2 and w’ 2 do not involve subtractions. To over- 
come this problem we have avoided a direct squaring of the radial velocities. In- 

2 

stead, we have used various interpolation formulae of v' r whose error is unbiased. 
Specifically, let i be an index for either a range or an elevation let A be the reso- 
lution and, let oti be some data We write, 
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(«.) 2 


a t [f^i + l + f a i-l — 6 a, + 2 ~ 6 a '- 2 ] + ^(A 3 ) > 

® i ^ [o, + i + tti-l ] + 0(A 2 ) , 

J a | ot , + 1 -t~ O ( A ) , 


a,«,-i + 0(A) , 


a t4-i a t-i + 0( A 2 ) , 


if r.h.s > 0 ; 
if r.h.s > 0 ; 
if r.h.s > 0 ; 
if r.h.s > 0 ; 
elsewhere 


(39) 


For convenience we define, 


2 i 2 1 1 

3 #t+l ' 3 — 1 g^i-f-2 g ^ t — 2 ? 

if (r.h.s)-a,- > 0 ; 

else 

a i 2 [ a '+! + Q «-l] i 

if (r.h.s)-®* > 0 ; 

else 

a H-l , 

if (r.h.s)-a* > 0 ; 

else 

<*i-l j 

if (r.h.s)-®, > 0 ; 

else 

[Oj+l^i-l] / , 

if other tests axe 

null 


( 40 ) 


Therefore, ( 39 ) is rewritten as, 


( a ,) 2 = a'"'a, + 0 (A n ), for n > 1 . ( 41 ) 

Recall now that, for the vertically pointing beam we have a spatially high 
resolution data in elevation (±2° from zenith with 0 . 2 ° increment). Therefore, for 
the vertical beam analysis the index z refers to elevation. For most parts the hori- 
zontal variance analysis is confined to elevations smaller than 20° and the index i 
then, refers to range. At any rate the error in estimating ( v'* ) is from ( 41 ) 
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s ( v 'r ) i = s [«);"'] «)i + « tw).-i «)i“' + < [(oj"] * i«)i] + o(A") . ( 42 ) 


It is worth noting that, 


(<(«' r ), v5(v')!”'> = 0 but ; (S(v' r )iS(v !0i) # 0 . (43) 

It is apparent from (39)-(43) and the fact that the errors should not be corre- 
lated with the “exact values” that, using the approximate form for the square of 
the radial velocity would result in an unbiased errors i.e., 


(*(«?)«) =0 + O(A>), for n > 1 


(44) 


The error covariance matrix of 6 (v 1 * ) can also be calculated. While the flux 
errors covariance matrix is diagonal this latter quantity is not, because, the in- 
terpolation induce correlations. However, the matrix is still quite sparse. Each 
observational error is correlated only with its two immediate neighbours. Further 
more, a straight forward calculation of the error correlation between two adjacent 
neighbours, and utilizing the various remarks about correlation (or lack of it) be- 
tween various quantities would reveal that these off-diagonal terms are smaller 
than the diagonal one by a factor of | cor [( . Where cor [ ] is the cor- 

relation coefficient (< 1) between two neighbouring observations. Therefore, we 
proceed with error analysis ignoring the effects of the off-diagonal terms. 

Regarding the standard deviations of 6 we utilize (40) — (44) and the 

supposed lack of correlation between errors and “exact’ 1 observations (25) to find 
that for an individual measurements, 
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a{8v’') = 2a(v' r )a(6v ' r ) . 


(45) 


And for N measurements, 


cr(8v' r ) = (2jN l ^<r{v' r ) cr{8v' r ) . (46) 

We can apply (46) directly to find (for vertically pointing beams) the stan- 
dard deviation of the error of the vertical velocity variance namely, 


a(8w’ 2 ) = (2/N l t^cr{v , r )<j(8v ' r ) ; for the vertical beam. (47) 

We can take a 2 (v' r ) = iu' J « 1 m 2 s -2 . As has already been indicated in the 
experimental design section for the vertically pointing beam we have, cr(8v ' r ) = 
±0.35 msi. For a single volume scan we have 40 observations which are only ±2° 
away from the zenith; we have 45 volume scans so that N = 1800 . Overall then 
cr(8w' 7 ) = ±0.016 m 2 s -2 ; a rather small number. 

Similar considerations can be used to estimate the errors in the third mo- 
ments of the vertical velocity with the result that, 


a(8w' 3 ) = ( 3/N 1/2 )w' 2 a{8v ' r ) . (48) 

Substituting appropriate values we get a(8w' 3 ) = 0.024 m 3 s -3 . 

We can also apply (46) to estimate the error in the horizontal wind variance. 
The following considerations are pertinent: 

First for heights larger than 600 m one can append the l.h.s of (10) and (11) by 
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the vertical velocity variance obtained from the vertically pointing beam. Second, 
for heights less than 600 m the horizontal variance analysis can be constrained 
so that observations with an elevation angle greater than 20° axe excluded. For 
a 600 m height this means that elevations in the interval 5.74° < © > 18° axe 
included. If we neglect the vertical velocity variance contributions in (10) or (11) 
the bias will be in the order of w' 2 sin 2 (^(18 + 5.74)) . For vertical velocity vari- 
ance in the order of 1 m 2 s -2 we get a small bias of the order of 4 cm/s. Third, 
by adding (30) to (29) we can eliminate the momentum flux from the equations. 
Taking the above into account we can follow steps similar to (31) — (37) and ob- 
tain, that 


a(6u' 2 ) — [2 1 / 2 N j / 2 / cos 2 (0)1 <t(u') a(8v' r ) . 


(49) 


Likewise, 


a(6v' 2 ) = |2 1/,2 iV 1/f2 / cos 2 (@)j<r(t;|.)o r (^?;|.) . (50) 

Typical values for N in the horizontal variance analysis are in the order of 
4000 for low heights to 2000 for z ~ 1.5km. We take 0 ~ 10 , a(6v' r ) = 0.7 ms -1 
and &(v' r ) — 2 ms -2 and obtain that, a ^(horizontal variance 0 0.05 m 2 s 2 

To test the effect of the interpolation formula (39) we have conducted exper- 
iments comparing the use of (39) vs. straight forward squaring. We find, that the 
momentum flux estimates are not sensitive to such variations. This is in line with 
the finding of the previous section where, we have shown that the flux estimate is 
not biased. By contrast employing (39) have resulted in variance reduction of the 
order of 1 — 2 m 2 s -2 . 
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6. Principal Features 

a. Mean properties 

Using our analysis technique we find that the mean surface wind is essen- 
tially from the south (direction ~ 185°). We have chosen a right-handed Carte- 
sian coordinate system such that, the positive x-axis is the direction of the mean 
surface wind. We call the x-component of the wind u and the y-component v. 

The dependence of u in height for the two case studies is shown in Fig. 3. The 
wind at the lowest 500 m appears to be well mixed; a southerly jet is prominent 
at ~ 1 km. The dependence of v in height is shown in Fig. 4. By construction v 
is zero at the surface. It is becoming more westerly with height up to an height 
of ~ 0.9 km for the earlier sounding and ~ 1.1 km for the later sounding, then 
becoming more easterly (at 1.3 km and 1.5 km respectively) where it turns more 
westerly again. All the above features are in good agreement with hourly releases 
of rawinsonde (Brutsaert, FIFE Document, ABL-1). In Fig. 5 we present poten- 
tial temperature soundings (Brutsaert, ibid.) starting at 15:45, 16:59 and 18:54 
UTC respectively (i.e., before, during and after our Lidax measurements). The 
growth of the mixed layer (i.e., the region where the lapse rate is close to adia- 
batic) during the two hour period is apparent. An exact definition of the “top” of 
the mixed layer is difficult. By one criteria we can place the “top” to be at 700 m 
for the first two soundings and, at 1.3 km for the last one. By another criteria 
one may place the top at 0.7 km, 1 km and 1.3 km respectively; more in line with 
traditional mixed-layer models of the type discussed by Tennekes (1973). 

The mean vertical wind is displayed in Fig. 6. One profile is derived using 
only the vertically pointing beam data. The second profile is estimated from data 
where the elevation angle is in the neighbourhood of 30 < O < 42. This latter pro- 
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file is used to define the mean radial velocities (9) and the fluctuations needed for 
the momentum flux calculations (32). While, the w estimate from the vertically 
pointing beams are less prone to observational errors, the statistics of the inclined 
beam may be more representative. It is reassuring to find essentially zero vertical 
velocities above the mixed layer. The unusually large mean vertical motion in the 
mixed layer is perplexing. It could be due to deficiencies of the Taylor hypothesis. 
Thus, long term sampling (1 hr) of the vertically pointing beam data may not ad- 
equately represent the eddies that are either too slow or too large. We have also 
attempted to obtain mean vertical motion by vertically integrating the continuity 
equation from the bottom up [(7-a), (7-b), (8)]. While the integration reveals neg- 
ative motions at the lowest 600 m there is still no satisfactory matching at 600 m 
between the vertical integration result and the vertically pointing beam data. 

b. Momentum fluxes 

The vertical profiles of u'w' are shown in Fig. 7 and, those of v'w' are dis- 
played in Fig 8. At the mixed layer the sign of the fluxes are in agreement with 
the eddy viscosity concept, i.e. 


u'w' = — K\dujdz ; and; v'w' = — Kidvjdz ; with; (7v i . A' 2 ) > 0 . 

The above is true provided that, the shear is interpreted not locally but, as some 
kind of vertical average from the surface layer to the top of the mixed layer. Tak- 
ing local values of the shear may lead to unrealistically large values of Aj and A 2 
when the wind does not change with height (Figs. 3 and 4). 

In the mixed layer the u'w' values are negative; near the surface this implies 


retardation of the surface wind. The values of v'w' on the other hand are posi- 
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tive; this implies retardation of the (negative) v component (Fig. 4). For compar- 
ison purposes we also display in Fig. 7 the aircraft estimate of u'w' at z = 50 m 
(McPherson et al., FIFE document, ABL-6). The Lidar and aircraft values are 
consistent even though the lowest level where Lidar derived fluxes are feasible is 
450 m. It should also be mentioned that with the exception (noted further be- 
low) the flux values are insensitive to the method by which w is calculated. The 
present results are obtained by using the w estimated from the inclined beams 
data (30° < 0 < 42°) data and then, subtracting it (and the mean horizon- 
tal wind) from the radial velocities. But, similar results are obtained if w is esti- 
mated from data where the beams are vertically pointing (±2°). The only excep- 
tion is that, for the 16:11-17:11 UTC period a local minima (not shown here) is 
indicated around z=800 in if w is estimated from vertically pointing beam data 
whereas, none seems to be apparent (Fig. 7) if the present w is used. 

The flux values of u'w' in the stable layer above the jet do not support the 
eddy viscosity concept; it continues to be negative even though the wind above 
the jet is decreasing with height. Furthermore, unusually large values are ob- 
served for both case studies (e.g., u'w' ~ —2.5 m 2 s -2 at z ~ —1.4 km for the 
17:29-18:20 UTC case study). To get further insight on the nature of the per- 
turbations we have examined time series of the vertically pointing beam data 
at various observational levels (600 m < z < 1950 m). The time resolution is 
36 s. Fig. 9 displays the time series for the period 17:29-18:20 UTC. Fig. 10 is 
the spectral analysis of Fig. 9. In the stable layer (say z > 800 m) the spectra 
has two distinct peaks; one with period T=54 min (the first Fourier mode) and 
the other with period of 10.8 min (the fifth Fourier mode). The significance of 
the first Fourier mode is not clear. According to standard time series analysis 
(e.g., Jenkins and Watts, 1968) the spectral amplitude of the first Fourier mode 
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may be influenced by artifacts such as the method of detrending, an inadequate 
sampling, the assumption of periodicity, etc. The peak of the fifth Fourier mode 
( 10.8 min) is deemed not to be susceptible to the above artifacts. It is interest- 
ing that the 10.8 minutes period correspond to an angular frequency (w = 2n/T) 
of the order of 10 -2 s -1 , not too far from the standard atmospheric value for the 
Brunt Vaisala frequency \{g / 6){d6 / dz)\ . A spectral peak at the Brunt Vaisala 

frequency has also been reported in numerous other studies of vertical velocity 
spectra of stable atmospheres (e.g., Ecklund et al., 1986 ; Gage, 1990 ). From the 
wind and temperature sounding we have also estimated the Richardson number 
R, defined viz., 


_ (jmde/dz) 

(du/dz) 2 ' 

In the stable layer we have found that, the Richardson number is in the range 
of 1 - 2 ; significantly above the critical Richardson number (= 1 / 4 ) below which 
turbulence is dominant. Furthermore, a Richardson number greater than unity is 
a necessary condition for waves absorption at critical layers. 

At present we interpret the large counter-gradient fluxes as a manifestation 
of the interactions between the convection in the mixed layer and waves in the 
stable layer. In particular, for reasons discussed further below we believe that 
the large fluxes are related to the existence of a so called “critical layer ’. Previ- 
ous studies (Kuetner et al., 1987 and Clark et al., 1986) have already identified 
the convection in the mixed layer as an important energy source for waves in the 
troposphere. To further proceed, let us examine the equation that describes the 
vertical structure of two dimensional linear gravity waves in a stable atmosphere 
is (e.g., Gossard and Hooke, 1975 p. 456) 
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W tz + 


TV 2 

(c — u ) 2 


^:z 

(c-u) 



W = 0. 


(52) 


Here c is the speed of the wave, W is the amplitude of the vertical wave velocity 
and k is the angular wave number (c = u/k ; with u> the angular frequency). Dif- 
ferentiation with respect to z is denoted by a subscript z and other symbols have 
been defined above. If a situation occurs where, the mean wind (u) at a layer of 
height z is equal to to the speed of the wave then, that layer is termed “critical 
layer”. Mathematically the differential equation (52) becomes singular. 

A detailed quantitative discussion of critical layers is beyond the scope of this 
paper but few remarks are in order. When a wave reaches a critical layer and the 
Richardson number is greater than 1, the wave energy is presumed to be com- 
pletely absorbed. Even if (due to say non-linear effects) reflection does occur its 
effect must be limited since, the wave can not significantly penetrate into the adi- 
abatic layer below. Overall then, if the predictions of linear theory are qualita- 
tively correct the presence of a critical layer should result in a build-up of wave 
energy and momentum in that layer. Of course, once the wave amplitude grows 
non-linear effects take place. As long as the flow can be characterised by a single 
frequency (e.g. the Brunt Vaisala frequency) the effect of the non-linearity is to 
steepen the wave front* and to diminish the gradients at the back. This should 
enhance the ‘‘downward” transport at the front and diminish the “upward” trans- 
port in the back so that on the average the “downward” transport is enhanced. 
Another consequence of non-linearity is" the breakdown of the so called Eliassen 
and Palm flux theorem which states that for linear gravity waves d/dz^u'w') = 0. 

* Here “front” and “back” are defined in a frame of reference moving with the 
mean wind; the terms “upward” and “downward” are used in the context of a 
right-handed coordinate system with the x-axis pointing in the direction of c — u. 
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Our data (Fig. 7) indicate large momentum flux divergence. 

To precisely establish the existence (or lack of) a critical layer one needs to 
independently measure the mean wind and the group velocity of the wave. This 
latter measurements are not feasible with the present scanning. Therefore, the 
above discussions must be regarded as speculative; however, if one accepts the 
proposition that the energy source is the PBL convection and, the Brunt VaisaT’a 
frequency is dominant, then it is possible to establish a bound for the phase ve- 
locities and see whether they are in a range which can overlap with that of the 
mean wind. From our spectra of the horizontal velocities (Figs 15 and 16) we 
have inferred that, the horizontal scales of the energy containing eddies in the 
PBL are in the range of at least 600 m - 6 km. This should excite waves in the 
stable layer with wave length A also in that range. Taking the Brunt Vaisala fre- 
quency N to be ~ 10 -2 s -1 we get (from c = u/k and k = 2 tt/X) that the phase 
velocities are in the range of 1 ms -1 — 10ms _1 . From Fig. 3 (the mean u wind) we 
establish that the velocities in the stable layer are in range of 8 ms -1 to 19 ms -1 
with an indication of further decrease with height. Thus, an existence of a criti- 
cal layer seems plausible. We suspect that such critical layers may be present in 
other cases where, the Richardson number is above unity and, the energy con- 
taining eddies in the mixed layer have lengths scales large enough, to excite waves 
whose phase velocities are comparable to the mean wind aloft. 

c. Variance 

The vertical profile of the variance of w w' 2 is shown in Fig. 11. We have 
used the vertically pointing beam for these calculations. Only the 17:29 - 
18:20 UTC is shown here because of some uncertainties in the pointing accu- 
racy of the vertically pointing beam for the earlier time. The variance gener- 
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ally decreases with height but a relative maxima near the critical layer is ob- 
served. Because vertically pointing data is not available at heights less than 
~ 0.52, ; ( Zi the height of the mixed layer), the well known maxima at ~ (0.2z, — 
0.4z;) is not detected. 

The vertical profiles of the variance of u ( u ' 2 ) is shown in Fig. 12 and that 
of of v' 2 is displayed in Fig. 13. For most part the profiles are consistent with 
the established conceptual picture of the PBL (e.g., Caughey, 1984). In the sur- 
face layer the variance is decreasing with height; in the mixed layer the horizon- 
tal variances seems to be well mixed as well. A notable departure from classi- 
cal results is in the vicinity of presumed critical levels. The horizontal variances 
seems to have local maximas when, local maximas of absolute values of momen- 
tum fluxes axe present. In the part of the stable layer above the critical layer, the 
fluxes diminish with height in agreement with other observations of the stable 
layer (e.g., Lenschow et al. 1989). The variance of u also diminishes with height 
at the top of the domain, but the v variance does not. This latter phenomenon is 
probably associated with the transient nature of the v wind in that region (Fig. 

4). Later soundings (not shown here) has revealed that the v wind (approxi- 
mately east-west in the present coordinates) to be consistently more westerly with 
height. 

d. The kinetic energy dissipation, and the surface heat flux. 

Recall from the discussions in section (4. c) the balance equation for the ver- 
tical velocity fluctuations in the form (18) rewritten here as (53) viz.. 


d . 1 — 3 , 
& ( 2 “ ) = 


1 ( < dP ' \ e - 9 


+ — w'O 1 


(53) 
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As has already been discussed the pressure term may be neglected near the 
surface therefore, the heat flux may be calculated as a residual balance between 
the the gradients of the third moments and the kinetic energy dissipation e. Re- 
call also that, the kinetic energy dissipation can be estimated from the spectra 
and (19) provided that an inertial «~ 5 / 3 law is established. 

Fig. 14 is a profile of the third moments of the vertical velocity w ' 3 . The 
positive values at heights 600 m < z < 800 m is encouraging. On the other hand 
the negative values up to an height of 1200 m axe perplexing. If it is real it could 
be due to the steepening of the waves or the or skewed negative motions in the 
entrainment zone. On the other hand it could be due to sampling; and the inade- 
quacy of relating turbulence parameters to the statistics of time series in a single 
spatial location. 

Figs. 15 and 16 are the longitudinal spectra of the u and v velocity correla- 
tion in the surface. All the beams in the lowest 100 m of the domain are included 
Each beam is spectrally analysed and than all the spectras in the entire time pe- 
riod (17:29 - 18:20) are combined. Prior to evaluating the Fourier coefficients the 
data is subjected to an high pass filter to remove wave length larger than 5.4 km. 
To minimize aliasing a low pass filter with weights (1/4, 1/2, 1/4) is also applied. 
This removes all wavelength shorter than 2A, where A is the range resolution. In 
both the u and v spectra an inertial -5/3 power law is clearly established. As the 
wave lengths get shorter, some attenuation due to the application of the low pass 
filter is apparent. The saturation at the high wave numbers is almost surely due 
to noise. 

From the amplitude of the power spectra at the inertial range we have es- 
timated that e — 5.02 x 10 -3 m 2 s -3 . By assuming that w' 3 (z = 0) = 0 
and, using the next two data points in Fig. 14 to fit a parabola we obtain 
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2 = 25 m) = 1.43 x 10 -3 m 2 s -3 . From this and using p 0 = 1.275 kg m -3 

, C p = 1004 J deg -1 kg -1 and Q 0 /g = 30degs 2 m -1 we get a surface heat flux 
value ~ 100 Wm" 2 . This is in apparent agreement with the average estimate of 
the surface stations. Nevertheless (as discussed above) the contribution from the 
gradient of the third moments is uncertain. 

6. Conclusions 

1. With respect to the analysis technique we have shown that for an horizon- 
tally homogeneous turbulence the mean wind, its variance and the vertical 
momentum fluxes can be deduced from a single Doppler Lidar data. While 
the accuracy of individual radial velocity is only ±0.7 ms -1 we have ex- 
ploited the large sample size to reduce the errors considerably. For the pur- 
pose of estimating the first moments the error reduction is in essence a sim- 
ple averaging procedure. In calculating higher moments non-linear effects 
usually modify the error field so that it becomes biased and correlated. Un- 
der these circumstances simple averaging does not reduce the error’s vari- 
ance. We have developed effective procedures to remove bias and reduce 
error correlation. Taking into account the large sample size we have shown 
that if the above mentioned procedures are utilized, we can obtain estimates 
of the variance and momentum fluxes with sensitivity to observational error 
which is at most ±0.10 m 2 s -2 . Regarding hardware improvements, the most 
serious problem is the stability of the reference frequency. Another serious 
limitation is the range resolution (=150 m). For PBL studies a more desir- 
able resolution is 25 m. New solid state Doppler Lidars could address both 
issues and are currently under developments in WPL and elsewhere. 

2. Spectral analysis of the horizontal velocities has in a natural way produced a 
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k - 5 / 3 power law at the inertial range. We view this as another indication of 
the viability of our analysis technique. From the power spectrum we have de- 
duced the kinetic energy dissipation. Using the statistical form of the Navier- 
Stokes equation the surface heat flux has been derived as a residual balance 
between the kinetic energy equation and the vertical gradient of the third 
moment of the vertical velocity. Our estimate of the surface heat flux is con- 
sistent with the average estimate of the surface stations. This latter result 
must be considered preliminary, because at present, the third moment of the 
vertical velocity cannot yet be reliably estimated at heights which are below 
600 m. 

3. The momentum fluxes in the mixed layer are in qualitative agreement with 
the eddy viscosity concept. However, in the stable layer the fluxes axe 
counter-gradient. The Brunt Vaisala frequency is the dominant one in the 
stable layer and the Richardson number is greater than one. We view this 
as an indication that the transport is accomplished by gravity waves. In the 
stable layer negatively large flux values ~ —1.5 m 2 s -2 are concentrated in 
thin layers where also, the flux divergence is large; in contradiction to the 
Eliassen-Palm flux theorem which, predicts (in the linear case and in the ab- 
sence of sources or sinks) that the flux divergence is null. Our interpreta- 
tion of the large flux values is that they are due to an absorption by a critical 
layer and the non-linear effect of steepening of gravity waves. 
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FIGURE CAPTIONS 


Fig. 1 RHI display in the direction of the mean surface wind. Note the southerly jet at 1 km. 
Refer to text for further explanations. 

Fig. 2 RHI display in the direction normal to the mean surface wind. Note the eddy structure in 
the mixed layer. Also, the wind becomes more westerly with height 

Fig. 3 Profiles of u with height for the two case studies. A right handed coordinate system is 
used such that the positive x-axis is the direction of the surface wind. The surface wind 
direction at 17:29-18:20 UTC period is -185° while the direction at the other period is 
-189°. 

* 

Fig. 4 Profiler 6f v with height for the two case studies. Note by construction v =0 at the 
surface. Refer to Fig. 3 and the text for further explanations. 

Fig. 5 Vertical profiles of 0 (the potential temperature) from nearby balloon measurements at 
15:45, 16:59 and 18:54 UTC, before, during and after our measurements period. 
Courtesy of W. Brutseart. 

Fig. 6 Vertical profiles of w with height for the period 17:29-18:20 UTC. One profile is 
obtained using the vertically pointing beams and the other from inclined beams. 

Fig. 7 u'w' profiles for the two case studies. Refer to Figs.3, 4 and the text for explanations of 
the coordinate system. The aircraft data is for 17:29-18:20 UTC (Courtesy of R. 
Desjardins and I. MacPherson). 

Fig. 8 v ' w ' profiles for the two case studies. Refer to Fig. 3 and the text for explanations of the 
coordinate system. 

Fig. 9 Time series of the vertically pointing beam Doppler velocities for various heights and for 
the period 17:29-18:20 UTC TJie horizontal axis is time in hour. The time resolution is_ 
36 s. The amplitude in ms* 1 is indicated on the vertical axis, while height is indicated by a 
series of horizontal lines. 


Fig. 10 Spectral analysis of Fig. 9. The periods T n (in minutes) are indicated in the horizontal 
axis. Also indicated is the Fourier mode n=T/T n with T the longest resolvable period 
(=54 min). 

Fig. 1 1 Vertical profile of the variance of w ( w' 2 ) with height for the period 17:29-18:20 UTC. 

Fig. 12 Vertical profile of the variance of u ( u' 2 ) with height for the periods 16:11-17:10 UTC 
and 17:29-18:20 UTC. Refer to Fig.3 and the text for explanations of the coordinate 
system. 

Fig. 13 Vertical profile of the variance of v ( v' 2 ) with height for the periods 16:1 1-17:10 UTC 
and 17:29-18:20 UTC Refer to Figs.3, 4 and the text for explanations of the coordinate 
system. 

9 

Fig. 14 Profile of the third moment of the vertical velocity w' 3 derived from the vertically 
pointing beams. The observation period is 17:29-18:20 UTC 

Fig. 15 Longitudinal spectra of the u velocity correlation at the surface. K=n2it/X ma3C ; 
n=l,2,...,N. 2^,* is the longest resolvable wavelength (=5.4 km) N=X max /A with 
A=150m, A is the spatial resolution. The -5/3 fit is also indicated. 

Fig. 16 Longitudinal spectra of the v velocity correlation at the surface. Refer to Fig. 15 for 
explanations. The -5/3 fit is also indicated. 
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